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We present a local control scheme to construct the external potential v that, for a given initial 
state, produces a prescribed time-dependent density in an interacting quantum many-body system. 
The numerical method is efficient and stable even for large and rapid density variations irrespective 
of the initial state and the interactions. The method can at the same time be used to answer 
fundamental ^-represent ability questions in density-functional theory. In particular, in the absence 
of interactions, it allows us to construct the exact time-dependent Kohn-Sham potential for arbitrary 
initial states. We illustrate the method in a correlated one-dimensional two-electron system with 
different interactions, initial states and densities. For a Kohn-Sham system with a correlated initial 
state we demonstrate the interplay between memory and initial-state dependence as well as the 
failure of any adiabatic approximation. 
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The v-representability question is one of the outstand- 
ing problems in density functional th eory (DFT) [J-Q. 
In time-dependent DFT (TDDFT) the question 

is whether, for a given initial state, there exists a lo- 
cal external potential v that yields a prescribed den- 
sity by solution of the time-dependent Schrodinger equa- 
tion (TDSE). In the case of a non-interacting system, 
f-representability amounts to the existence of a Kohn- 
Sham (KS) system [l], 0|. The KS system has played a 
major role in the study of correlated many-body systems 
as it allows for the treatment of interacting systems in an 
effective one-particle framework. This feature greatly re- 
duces computational costs [12[ and (TD)DFT has hence 
been one of the leading methods in electronic structure 
theory [l3|. In practice the accuracy of the method is 
limited by the approximate nature of the density func- 
tional that are used. In TDDFT the most commonly 
used density functionals are based on the adiabatic ap- 
proximation in which the KS potential only depends on 
the instantaneous density. These functionals can, how- 
ever, fail in important cases [9| and therefore there is a 
great need for better functionals. To develop and bench- 
mark such new functionals the availability of exact time- 
dependent KS potentials is highly desirable. Although 
such potentials can be constructed in special cases [14jJl5j 



no general practical scheme has been available so far. In 
this Letter we provide such a scheme based on a rece ntly 
introduced fixed-point formulation of TDDFT [Io|, flq . 
It is at the same time an efficient local control scheme 
based on the densit y w hich augments other important 
control methods [i7H2q| already of extensive use in laser 
physics, quantum optics |2l| and the physics of ultra- 
cold gase s [22I. T he scheme is closely related to existing 
methods [19|, [20| but targets a spatially extended quan- 
tity instead. It is applicable to general interactions and 
initial states and can deal with fast and large density 



changes. We demonstrate the approach for an interact- 
ing one-dimensional two-electron system with different 
interactions, initial states and densities. For a KS system 
with a non-separable initial state we illustrate the con- 
nection between memory and initial state dependence as 
well as the failure of any adiabatic approximation. 

The global fixed point method. We consider a TV- 
electron system with a time-dependent Hamiltonian 
H(t) = f + V(t) + W, where f is the kinetic energy, V(t) 
the time-dependent external potential and W the many- 
body interaction (which may even be time-dependent). 
The expectation values n(rt) and j(rt) of the density and 
current operators (atomic units are used throughout) 

TV 
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satisfy equations of motion given by 

d t n(rt) = -^-j(rt) 
d t j(rt) = -n(rt)^v(rt) + Q(rt). 

Here the internal local force Q(rt) is defined by 
Q(rt) = -i<*(M,t)|[j(r),f + W]\V([v],t)) 



(1) 
(2) 



where |^([v],t)) is the time-dependent many-body state 
obtained from the TDSE with potential v and given ini- 
tial state. Eqs.(pQ) and (j2j) imply 

- ^ • (n(rt)^v(rt^j = q([v],rt) - d?n(rt). 

where q([v],rt) = — ^ • Q(rt) is regarded as a functional 
of v through the state \^([v],t)). For a fixed density and 
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initial state this is an implicit equation for the potential. 
To solve this implicit equation we define an iterative se- 
quence v k of potentials by the iterative solution of 

- ^ • (n(rt)^v fc+ i(rt)) = q([v k ],rt) - d 2 n(vt). (3) 

In previous works [Io|, Esj we proved, for general initial 
states and interactions, that under mild restrictions on 
the density the sequence v k converges in Banach norm 
sense to a potential v which is both fixed-point of the 
equation and produces the prescribed density n. 
Although the fixed point method itself is well-defined it 
is highly non-trivial to develop a stable numerical algo- 
rithm. To do this we found it advantageous to make 
explicit use of also the current (still being a functional of 
the density). This is most easily done for one-dimensional 
systems since the continuity Eq. © can be integrated an- 
alytically. We find j(xt) = j(at) + j(xt) where 



f 

J a 



dyd t n(yt), 



(4) 



and where a is an arbitrary point. For this reason, and 
for simplicity of presentation, we restrict ourselves to the 
one-dimensional case in this Letter. We first show how 
we can eliminate the quantity q from our equations. By 
integrating Eq.© and using Eq.© we obtain 

- n(xt)d x Vk+i(xi) =d t j(xt) - Q([v k ], xt) + ci(t), 

where c\(t) is an integration constant. From Eq.© for a 
system with potential v k we then find 

- n(xi)d x v k +i (xt) = d t [j(xt) - j([v k ],xt)] 

- n([v k ],xt)d x v k (xt) + ci(i).(5) 

To obtain an equation that is only dependent on densities 
we can use Eq.© and Eq.© to find 

- n(xi)d x v k +i(xi) = / dyd 2 [n([v k ],yt) - n(yt)] 

J a 

- n{[v k ],xi)d x v k {xt) + c 2 {t), (6) 

where C2(t) is a new constant. While mathematically 
equivalent Eqs.© and Eq.© are not numerically equiv- 
alent as their discretizations on a space-time grid gen- 
erally differ. In practice, it is therefore advantageous to 
use 

- n(xt)d x v k+1 (xt) + n{[v k ],xt)d x v k (xt) = 
(l-/i) / dyd 2 [n([v k ],yt)-n(yt)] 

J a 

+lidt [j(xt)-j([v k ],xt)] +c(t) (7) 

as follows immediately by multiplying Eq.© by /i and 
Eq.© by 1 - /i and adding the results. Here fi is a 
parameter at our disposal and c(t) is a new constant. 



This equation defines an iterative procedure to deter- 
mine v k+ i from v k . The constant c(t) in this equation is 
uniquely determined by the spatial boundary conditions 
on v k +\ (and hence depends on k). When v k — >• v then 
n[v k ] — >• n[v] and j[v k ] — )• j[v] and Eq.([7j) implies that 
c(t) jidtj(at). Since we also obtain I^Qu], t)) after 
convergence we can calculate any observable, and in par- 
ticular the current j(xt). This is an explicit realization 
of the Runge- Gross result [23j that any observable is a 
functional of the density and the initial state. 
Numerical procedure. The iterative method based on 
Eq. ([7]) should be implemented stepwise in time for 
high efficiency. We use a midpoint based time-stepping 
method which uses the midpoint potentials v{xt n ) = 
v(x, \ {t n -\ + t n )) to propagate the wave function on a 
time-grid with time-points t n . For this we implemented 
the Split Operator and Lanczos method . Let us now 
suppose that we have obtained v(xt m ), and hence the 
|^([v], tm)) giving the required density n(xt m ) for m < n. 
Then to determine v(xt n +i), and hence |^([f], t n +i))> we 
define an iterative procedure in which we guess an initial 
potential and loop over potentials v k (xt n +i) until we con- 
verge to the desired v(xt n +i): 

1. Use v k (xt n +i) to calculate |^([^], t n +i)) fr° m 
\V([v],t n )) by time-stepping. 

2. From ^([^J, t n +i)) calculate n([v k ],xt n +i) and 
j([v k ],xt n +i). 

3. Calculate v k +i(xt n +i) from 

-n(xt n+ i)d x [v k +i(xt n +i) - v k (xt n+1 )} At 2 



dy [n([v k ],yt n +i) - n{yt n+1 j\ 
■ j([v k ],xt n +i 



(8) 



+ BAt[j(xt n+1 ) 

where At = t n+ i — t n . 

Eq.© is obtained from Eq.© by a discretization w.r.t. 
time using only times t m with m < n + 1 for the deriva- 
tives. We further used the fact that we have already 
converged up to time t n and replaced n[v k ] by n on the 
left hand side of the equation since we found that this 
does not affect the convergence. The constants A and B 
depend on the discretization scheme and the fi of Eq.©, 
which effectively leaves the choice of their values at our 
disposal. The constant c in Eq.© depends on the bound- 
ary conditions and hence the geometry of the system. 
Below we will present examples for a periodic system. In 
that case the constant is determined by the periodicity 
condition v k (at) = v k (bt) on the spatial interval [a, b]. 
This yields the condition 



dx 



Af*dy [n([v k ],yt n +i) - n(yt n+1 j\ 



n(xt n +i) 



BAt [j(xt n+1 ) - j([v k ],xt n +i)] 
fi(xt n+ i) 



0. (9) 
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FIG. 1: (color online) The potentials that produce the pre- 
scribed densities m and ri2 (insets). Panel (a) m, A = 0, (b) 
n2, A = 0, (c) ni,A = 1, (d) ri2,A = 1. Note that in (a) we 
plotted minus the potential for better visibility. 



During time-propagation numerical errors can build up 
in the modulus and phase of the wave function. The 
density is determined by the modulus only while the cur- 
rent also depends on the phase. This implies, for ex- 
ample, that errors in the phase can lead to inaccurate 
currents while still producing an almost correct density. 
This tends to happen in the case that B = since in that 
case the procedure enforces the correct density without 
constraints on the current as can be seen from Eq.(j8]). 
The opposite happens in the case that A = 0. By taking 
nonzero values for A and B we control the accuracy of 
both the density and the current and hence the modulus 
and phase of the wave function. This suffices to stabilize 
the algorithm in most cases. Perfect stability is obtained 
by spatially smoothening the potentials Vk+i as they are 
obtained. Without smoothening the best algorithm is 
obtained when A dominates B for nonzero B while these 
values are not so important with smoothening (we used 
A = 1 and B = 0.5). We find that 5 iterations generally 
suffice to converge and that the precision of the potential 
is limited mainly by the time-stepping method for the 
wave- function (assuming a sufficient spatial resolution). 
By increasing the precision thereof almost arbitrary pre- 
cision can be achieved even when the density changes by 
orders of magnitude. 

Translating and splitting a given density. To illus- 
trate the algorithm we consider two electrons on a quan- 
tum ring of length L = 10 over a time period of length 
T = 20. We start by calculating the singlet ground state 
|^o) (which has a spatially symmetric wave function) of a 
(properly periodic) Hamiltonian with external potential 



and interaction w given by 
(x) = — cos 
w x 2 ) =Acos 



(2ttx\ 

27T (Xi - X 2 ) 



where A is the interaction strength. The ground state 
density is denoted by no(x). We then construct the (spa- 
tially periodic) time-dependent densities n\ and n 2 by: 



ni(xt) 
n 2 (xt) 

r(f) 



n (x - r(t)), 

- [n (x - r(t)) + n (x + r(t))] , 

L [ fnt 

— 1 — cos — 
2 V T 



The density n\ describes a situation where the initial 
density no is rigidly translated around the ring exactly 
once whereas the density n 2 describes a situation where 
the initial density no is split in equal halves ^no that 
are rigidly translated in opposite directions to rejoin at 
times \T and T. We have used our algorithm to calculate 
the potentials that produce these prescribed densities ri\ 
and n 2 via time-propagation of the initial state |^o) by 
the TDSE. This was done for the interaction strengths 
A = and A = 1. In Fig. [I] we present the corresponding 
potentials and densities (insets). We see large differences 
in the potentials for the interacting case (panels (c) and 
(d)) as compared to the non-interacting case (panels (a) 
and (b)). The convergence of our algorithm shows that 
the prescribed densities are indeed v-representable and 
that the algorithm can be used for density changes of 
orders of magnitude. 
Exact KS potential for a non-separable initial state. As 
a second example we construct an exact KS system, i.e. a 
non-interacting system having the same time-dependent 
density as that of an interacting reference system. For 
the KS system we also need to specify an initial state 
with the correct initial density no- This state does not 
need to be the KS ground state (the ground state of a 
non-interacting system with density no) as the Runge- 
Gross theorem [23[ allows for general initial states (see for 
further discussion [HI). Here we take the KS initial state 
to be identical to the true correlated ground state |^o) 
of the interacting system. As the interacting reference 
system we consider a system forever kept in the ground 
state |^o) of the previous example for A = 1. The density 
is therefore stationary and equal to no- Since |^o) is not 
an eigenstate of a noninteracting system the KS state and 
potential will in general be time-dependent, but in such 
a way that they still produce the static density no- We 
denote the KS potential by v s and the KS Hamiltonian 
is thus given by 



H s {i) 



1 



(<9? + <9f) +v 8 {x 1 t)+v a {x 2 t). 




FIG. 2: (color online) The explicitly time-dependent KS po- 
tential that keeps the density (inset) static for the correlated 
initial state |^o)- We stress that the potential is not periodic 
in time. 
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states where they are confined to the same region in space 
(t = 8), although the corresponding density is completely 
static. The wave function at these times as well as the in- 
termediate times t = 4 and t = 12 are in correspondence 
with the extreme values of the potential in Fig|2] We 
also see the exact potential cannot be an adiabatic func- 
tional of the density, and hence must have memory, as an 
adiabatic functional produces a static potential when we 
insert the exact density, in conflict with FigJ21 This can 
be illustrated further by choosing as (the spatial part of) 
the initial KS state a separable state of the form 



^o(x 1 ,X 2 ) = (t>(xi)(j){x2) 



(10) 



with 4>{x) = ^/tiq(x)/2. In this case the KS-potential v s 
is static and given by 



v s (x) 



1 dl^np(x) 

2 \/n Q (x) 



(11) 



up to an arbitrary constant. In this case the exact KS po- 
tential is static as would also have been predicted by any 
adiabatic approximation. This explicitly demonstrates 
the interplay between memory and initial states (26i l27|. 
Outlook. We presented a stable and fast algorithm to 
construct the external potential that, for a given initial 
state, produces a prescribed time-dependent density in 
an interacting many-body system. The method will be 
valuable for further development of density functionals 
and local control theory. Especially exciting is the possi- 
bility to use more advanced (multi-configurational) initial 
states in DFT in combination with existing and new ap- 
proximate functionals and to test them using our bench- 
marking algorithm. This can open up new possibilities 
for the study of strongly correlated systems within a DFT 
framework. 
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FIG. 3: (color online) 4 snapshots of |^ s (xi, X2, £)| 2 at times 
where the KS potential in Fig. ([2]) is extreme. Note that the 
electrons are well-separated at time t — but are confined to 
the same region at time t — 8. 



We have determined the time-dependent potential v s 
with our algorithm and displayed it in Fig|2j The square 
\$ s (x\, X2 : t)\ 2 of the corresponding KS wave function is 
displayed in Fig|3]at four times corresponding to extreme 
values of the KS potential. We see strong internal mo- 
tions in the wave function as it passes through states 
in which the electrons are well-separated (t = 0) and 
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